Search for periodic gravitational radiation with the ALLEGRO gravitational wave 
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We describe the search for a continuous signal of gravitational radiation from a rotating neutron 
star in the data taken by the ALLEGRO gravitational wave detector in early 1994. Since ALLEGRO 
is sensitive at frequencies near 1 kHz, only neutron stars with spin periods near 2 ms are potential 
sources. There are no known sources of this type for ALLEGRO, so we directed the search towards 
both the galactic center and the globular cluster 47 Tucanae. The analysis puts a constraint of 
roughly 8 x 10"^'* at frequencies near 1 kHz on the gravitational strain emitted from pulsar spin- 
down in either 47 Tucanae or the galactic center. 
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I. INTRODUCTION 



The majority of experimental searches for gravitational radiation have focused on the detection of burst signals, 
such as those emitted from the collapse of a massive star. There are compelling arguments that nearby millisecond 
pulsars can provide a detectable source of continuous (CW) gravity waves We will use the term "pulsar" in 

the following text to refer to a rotating neutron star, regardless of whether it is emitting detectable electromagnetic 
^ ■ radiation. 

It is well known that a symmetric object rotating about its symmetry axis does not emit gravitational radiation. 
Therefore, if a pulsar is to radiate gravity waves, it must either be asymmetric or precessing. For the purposes of this 
search, we have considered only asymmetric pulsars. The asymmetry can be rotationally induced or due to a "star 
quake" deforming the neutron star crust or due to some other process. By any mechanism, the amount of ellipticity 
produced in the neutron star is expected to be small, from a maximum of 10~^, which, for some pulsar models is 
[ the maximum supportable asymmetry of the neutron star crust ||] to less than 10^^. The latter number is obtained 
' by assuming the measured spin-down of certain millisecond pulsars is due entirely to the emission of gravitational 
radiation. 

I In this article we report the results of our first search for a CW signal in data taken by the ALLEGRO resonant 

^jpj gravitational wave detector at Louisiana State University during the first three months of 1994 Another search 
• • . has been performed by the EXPLORER detector team from the University of Rome, which has a similar detector 
. ^ ' s-^d is pursuing a different type of analysis |^ . 

, If there were known pulsars with the right spin period, they would be the natural target for any CW search. But 
5^ • there are only a small number of known pulsars with spin periods small enough (2.2 milliseconds) to match our 
[ antenna's 900 Hz frequency and none of those listed in the 1995 Taylor pulsar catalog has a spin period that 
matches our two narrow reception bands. We would ideally like to make an all sky search, but this is computationally 
very intensive. Given limited computing resources it nearly impossible to exhaust all of the obvious possibilities. 
All-sky search strategies have been proposed elsewhere . 

Therefore we have chosen another strategy and directed our first search towards two regions where one might 
reasonably expect to find a high density of CW sources. The first is the globular cluster 47 Tucanae (Tuc), and the 
second is the center of the galaxy. The celestial coordinates for both are listed in Table ||. We chose 47 Tuc because 
it is relatively close, and because of the large number of fast millisecond pulsars known to be concentrated in this 
globular cluster which suggests the possibility of more undetected neutron stars with short spin periods. 

The analysis assumes the simplest possible CW source, one that has a stable emission frequency in its rest frame. 
In other words it is not spinning down and has no companions to cause orbital doppler shifts. In this case only 
knowledge of the earth's motion and the antenna's reception pattern are needed to determine the modulation of the 
frequency and amplitude of the received signal. These modulations are highly dependent on the source location. 
At the frequency resolution of our search, if the actual source location is offset from the supposed location by as 
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little as 1.00' in right ascension and 0.25' in declination, the gravitational wave frequency will be mis-identified by 1 
bin 1^. However, the range of possible locations for detection of the signal (albeit with the frequency mis-identified) 
is significantly larger |10|. 



II. ALLEGRO 

The ALLEGRO gravitational wave detector |^ is located in the Physics Building at Louisiana State University in 
Baton Rouge, Louisiana (30°. 413 N lat, 91°. 179 W long). ALLEGRO consists of a resonant aluminum bar equipped 
with a resonant inductive transducer and a dc SQUID amplifier all cooled to 4.2 K. The bar is 60 cm in diameter and 
300 cm in length, with a physical mass of 2296 kg. The bar is oriented perpendicular to the plane of the great circle 
on the Earth that passes through Geneva, the location of the Rome Explorer antenna, and midway between Baton 
Rouge, LA and Stanford, CA. This orientation results in the axis of ALLEGRO being directed along a line 40°. 4 west 
of North. 

Vibrations of our resonant antenna produce a voltage out of the SQUID electronics which is proportional to the 
relative displacement of the bar and transducer. The coupled system of bar and transducer has two normal modes 
at roughly 896.8 Hz and 920.3 Hz. These are referred to as the "minus" and "plus" resonant modes respectively. 
ALLEGRO is most sensitive in a small bandwidth around each of these modes. The voltage out of the SQUID 
electronics is sent to a single lock-in detector which demodulates and low pass filters the signal. The reference 
frequency of the lock- in is set halfway between the normal mode frequencies of the antenna, thus shifting the frequency 
of the signal from near 1 kHz to near 10 Hz. The output consists of in-phase (x) and quadrature (y) components 
which contain the full amplitude and phase information. The demodulated data is then sampled 125 times a second 
and written to disk. 



III. TARGET SOURCE 



We describe the gravity wave in a reference frame which has its e^' axis aligned with the direction of propagation of 
the wave. This frame is referred to as the "wave frame" and is denoted by primed coordinates {x ,y ,z ). The most 
general form of the gravity wave in the wave frame is 



h+{t) h^{t) 

hyXt ) -h+{t ) 





(1) 



where /i+(t ) and hx{t ) are the amplitudes of the two allowed states of linear polarization, referred to as the plus 
and cross amplitudes respectively. 

In General Relativity one can compute the polarization amplitudes for a rotating, asymmetric neutron star using 
the quadrupole approximation ]ll[| . We choose for our model a neutron star with the three principle moments of 
inertia about the three principle axes fixed in the body frame, denoted by Ii, I2 and 13. The rotation axis is chosen to 
be along I3 and the neutron star is assumed to be deformed such that Ii ^ I2 ■ The ellipticity of the pulsar is then 
defined to be 



(2) 



The resulting gravitational radiation is emitted at twice the pulsar rotation frequency and the two polarization 
amplitudes are given by 



/i+(i') = hc{l + cos^ i) cos(27r/ot') 
hx(t ) = 2 he cos i sin(27r/ot ) 



(3) 



where the angle i is measured between the pulsar spin axis and the line of sight to the detector, is the "characteristic 
amplitude" [y of the incident wave, and fo = 2frot is the frequency of the gravitational wave, equal to twice the 
rotation frequency of the pulsar. This is not the most general pulsar model one could construct. If the neutron 
star deformation is not along a principle axis, then emission occurs at the rotation frequency and twice the rotation 
frequency [|l2| . If the star is precessing, emission occurs at the rotation frequency, and the rotation frequency plus or 
minus the precession frequency |l^] . For this work, we consider only emission at twice the rotation frequency. 
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The characteristic amphtud^ is given by 

K = -^eh {nfo? (4) 
c^r 

where r is the distance to the source, G is the gravitational constant, and c the speed of hght. Adopting a value of 
I3 = 10'^^ gcm^ (a reasonable estimate for a IAMq pulsar with a radius of 10 km), this can be written as 

he = 5.2 X 10--(-^) (-^) (^) (-A-r . (5) 
10*'' g cm^ 10 ° r 1 kHz 

For galactic sources (< 10 kpc), and assuming the maximum allowed pulsar ellipticity (e ~ 10^^) we see that the 
characteristic amplitude is of order 10^^^. As will be shown, this is of the same level as the measurements. 

IV. PHASE CONSIDERATIONS 

Reliable detection of a periodic gravitational wave signal depends on tracking the phase of the signal through many 
cycles. As given in Eq. ^ the polarization amplitudes of the gravity wave are pure sinusoids at an unknown frequency. 
Detection of such a signal would involve a single Fourier transform of the data, the resulting spectrum then being 
scanned for anomalous peaks. However, the phase of the gravity wave as observed at the detector at a particular 
observation time, relative to the emitted phase, depends on a number of factors: 1) intrinsic pulsar spin-down, 2) 
motion of the Earth within the solar system, 3) if it is part of a binary system, orbital motion of the pulsar, and 4) 
proper motion of the pulsar. This results in the initially narrow-band signal being spread out over many frequency 
bins, decreasing the signal to noise ratio of the single Fourier transform. As the expected signals are already weak, 
this seriously compromises one's detection capabilities. For the purposes of this work, we assume that the pulsar is 
solitary with negligible spin-down, and has no other accelerations with respect to the solar system barycenter (SSB). 
We also assume the proper motion of the pulsar is small such that it does not move significantly during the observation 
time. 

It should be noted that the spindown due to the emission of gravitational radiation at the level of sensitivity of this 
experiment would be quite significant. The lack of an actual target source has led us to choose the simplest case for 
the purposes of demonstration of the potential sensitvity of a real detector. With these assumptions, the phase of the 
carrier frequency is modulated only by motion of the Earth in the solar system, with the modulation being given by 
(see Fig. I) 

= ujo{ -r(<) • n + 1.658'"" sin g{t) ) (6) 
c 

where r{t) is the vector from the SSB to the center of mass of the detector at observation time n is a unit vector 
from the SSB directed towards the pulsar, c is the speed of light, and loq — 27r/o. The second term is due to the 
combined effect of time dilation and gravitational redshift due to the solar system bodies. It is periodic over a year 
with maximum delay of 1.7 ms. We note that the angle g{t) varies slightly from year to year. In 1994 its value was 

g{t)/2Tr = 356.60 + 0.98560 x D 

where D is the day of the year. Including this information into Eq. ^ we write the two polarizations of the gravity 
wave as (expressed in a frame parallel to the wave frame) 

h+{t) ^hc{l + cos^ i) cos(cjot + *(t) + *o) (7) 
hxit) = he (2 cos i) sin(wo^ + ^(t) + ^0) 

where $0 is an unknown, constant offset. 

Astronomical locations are usually given as right ascention and declination coordinates, denoted a, S respectively. 
These are defined in the the "celestial" coordinate frame (CC frame with coordinates {X, Y, Z)). This frame is centered 



'^New's characteristic strain is greater by a factor of 2 as they assumed that the rotation axis of the pulsar was along the line 
of sight to the Earth. 
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at the solar system barycenter (SSB) with ez along the Earth's rotation axis, ex and ey are in the Earth's equatorial 
plane with ex directed towards the intersection of the equatorial plane with the Earth's orbital plane (the ecliptic) 
at the vernal equinox. Right ascension is measured in hours of angle (12 hrs = tt radians) from the vernal equinox 
eastward along the celestial equator to the celestial object and declination is measured in degrees north (+) or south 
(-) of the equatorial plane. 

To calculate the dot product in the phase modulation term, we also expressed the detector position in the CC 
frame. This was done in two steps. First, the position of the center of mass of the Earth relative to the SSB was 
obtained using a commercially available software package from the U.S. Naval Observatory (MICA^). MICA uses the 
Jet Propulsion Laboratory DE200/LE200 ephemeris. It reports positions of the Earth's center of mass in a cartesian 
coordinate frame centered at the SSB. Coordinates are given to the nearest 10~^ astronomical units (AU), which is 
of order 10^ meters. This corresponds to 1/1000 of a wavelength for the signals of interest here, enabling an accurate 
tracking of the phase of the gravitational wave. 

Second, a GPS receiver was used to obtain the latitude and longitude coordinates of ALLEGRO, again with sufficient 
accuracy to track the signal. These were converted to cartesian coordinates in a frame centered at the Earth's center 
of mass and vectorally added to the center of mass positions to provide ALLEGRO'S position with respect to the SSB. 
The cartesian coordinates for ALLEGRO were converted to the spherical coordinates (r, t^, (5yi): r being the distance 
from the SSB to the detector, tg the local sidereal time, and 8a the declination of the detector. Using the spherical 
trigonometric formula for the cosine of the angle between two vectors, we write the phase delay as 

$(i) = 27r/o ^ [sin sin b + cos 8a cos 8 cosj{t) ] + 1.658'"'' sin g{t) ^ (8) 

where j = ~ a is the local hour angle. 



V. SIGNAL CONSIDERATIONS 



A passing gravity wave provides the largest fractional change in the length of a bar, and therefore the largest signal 
in a bar detector, when its direction of propagation is perpendicular to the bar axis and one of the polarizations of the 
gravity wave lies along the bar axis. In general, the gravity wave is incident to the bar with some arbitrary angle and 
polarization. We define the polarization angle as the angle between the bar axis (west of North) and the direction of 
the plus polarization of the gravity wave. 

We describe the detector in a coordinate frame whose origin is at the center of mass of the antenna, with the 
axis aligned with the local vertical and the e^ axis aligned with the bar axis. This frame is referred to as the "lab 
frame" and is denoted denoted by unprimed coordinates {x,y,z). 

For a gravity wave incident with an arbitrary orientation to the bar axis, it is only the component of the strain 
tensor along the bar axis which produces a detectable driving force on the bar. This force is commonly written as 

T{t) ^ ^flLe'Lxit) (9) 

where hxx{i) is the second time derivative of the strain component along the bar axis. The quantities /i and Le are 
the effective mass and length of the first longitudinal eigenmode of the bar, obtained by solving the elastic equations 
of motion for a long, thin cylinder. The effective mass is equal to half the physical mass of the bar (/i = 1148 kg) and 
Lf, = (4/7r^)L where L=3 m is the actual length of the bar. The force acting on the bar dominates the output as the 
force on the transducer itself produces a much smaller motion. 

To calculate the component of the gravitational wave, given by Eq. |l| in the wave frame, along the bar axis requires 
knowledge of the rotation matrix between the wave frame and the lab frame. As both source direction and detector 
locations are known, the rotation matrix is completely specified (up to an angle related to the unknown polarization 
of the gravity wave). We calculate the rotation matrix by first rotating the wave frame to the CC frame, and then 
the CC frame to the lab frame. The full rotation matrix from the wave frame to the lab frame is given by 
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4 



The CC frame is related to the wave frame by the angles {4> = a — Tr/2,9 = S + 7r/2,'0 = ipo), where we have 
used the Euler x-convention to define the axes of rotation. Using Eq. 4-46 of Goldstein , with the stated angular 
substitutions, the rotation matrix from the wave frame to the CC frame is given by 

((sin a cos tpo — cos a sin 5 sin -^o) (~ sin ipQ sin a — cos i/jq cos a sin S) (— cos a cos S) 
(— cos ^0 cos a — sin ipQ sin a sin S) (sin -00 cos a — cos ipQ sin a sin S) (— sin a cos 6) 
(sin ■00 cos (5) (cos -00 cos (5) (—sin (5) 

Again using the x-convention to define the rotation axes, the CC frame is related to the laboratory frame by the 
angles {(f> = ts+TT/2,9 = Tr/2 — l,ip = TT/2 + ri). is is the detector local sidereal time as before and I is the detector 
latitude. The particular value for the final rotation angle comes from the choice to define the axis pointing west of 
North along the bar axis. The rotation matrix from CC frame to lab frame is then 

((sin 77 sin — sinZ costs cos 77) (— sin 77 costs — sinZ sin ts cos rj) (cos rj cos I) 
(cos rj sin tg + sin I cos tg sin 77) (— cos 77 cos tg + sin I sin ts sin 77) (— sin 77 cos I) 
(cos Z costs) (cos/ sin ts) (sinZ) 

Using Eq. |l^, the full rotation matrix may be computed and h^x is the (11) component of 

h/ab — "^wave — dab^wave'^njave — >lab ' (-^-^) 

Considering only the component of the gravity wave along the bar axis, 

hbar{t) = (Rli - RI2) h+{t) + 2 Rn R12 (t) 
where the relevant components of the rotation matrix are 

Rii = cos tpo (sin 77 cos 7 -|- sin / cos 77 sin 7) 

+ sin tpoi— sin 77 sin 5 sin 7 -I- sin I cos rj sin S cos 7 + cos 77 cos I cos S) 

R12 — cos ipo{~ sin 77 sin S sin 7 -|- sin I cos 77 sin S cos 7 + cos 77 cos I cos S) 
— sin t/jq (sin 77 cos 7 + sin I cos 77 sin 7) . 



Defining 



and 



we have 



and 



/+ (t, ri,j,S,l) = (sin 77 cos 7(t) + sin I cos rj sin 7(t))^ (12) 
— ( — sin 77 sin S sin7(t) + sin / cos 77 sin S cosj{t) + cos 77 cos I cos 5)^ 



fx (t, ri,^,6,l) = 2 (sin77COS7(t) + sin / cos 77 sin7(t)) (13) 
X (— sin rj sin S sin 7(t) + sin I cos 77 sin S cos 7(t) + cos rj cos / cos S) 



Rli - RI2 = cos 27/'o /+(t) + sin 27/.0 /x (i) 
2i?iii?i2 =cos20o /x(t) -sin27/'o/+(t) 



hbar{t) = [cos 200 /+ (t) + sin 27/^0 /x(t)]/7+(t) + [cos 200 /x(t) -sin 27/'o/+(t)]/7x(t) . (14) 
Using Eq. Q and Eq. |lj, the waveform is given by 

hbar{t) = Kil + cos^i) [cos 2^-0 /+ (t) + sin 200 /x(t)] cos( wo* + ^(0 + *o ) (15) 
-I- he (2cos7) [cos 20^0 /x(0 - sin 200 /+(t)] sin(u;ot + *(t) + <i>o ) ■ 

The reception patterns f+it) and /x (t) are shown in Fig. ^ for a signal from 47 Tuc. Writing the time-dependent sine 
and cosine terms in Eq. |l^ as exponentials and defining 

F+(t)=cxp[j$(t)]/+(t) (16) 
Fx(t) = cxp[j$(t)]/x(t) 
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the expression for the gravity wave strain along the bar axis is 

hbar{t) = (1 + cos^ i) exp(j$o) [cos2tPoF+{t) + sin2V'oi^x (t)] exp{jujQt) (17) 

+ ^hc{l + cos^i) exp(-j$o) [cos2i/'o^*(i) +sin2iAo^'xW] exp(-jwoi) 
j 

- -hc{2 cos i) exp(j$o) [cos 2?Ao ^x (t) - sin 2?Ao ^+ (^)] exp(jLJoi) 
j 

+ -hc{2 cos i) exp(-j$o) [cos 2'(/'o -F* (t) - sin 2ipo (t)] exp(- jcjoi) • 

We are now in position to calculate the anticipated signal as it would appear in the ALLEGRO data stream. This is 
most clearly presented in the frequency domain. Fourier transforming Eq. ^, the force produced on the bar is 

^ -^^iLeUj'^ hbariul) ■ (18) 

This driving force produces motion of the transducer 

H{oj) = G{uj) J-{U!) = ——IJLLi, L0^G{L0)hljar{^) 

where G{uj) is the transfer function which relates transducer motion (or equivalently voltage) to the applied force [ p7[ . 
We note that the overall calibration is contained in the transfer function. Defining a new transfer function (strain to 
transducer motion) 

Gf{L0) = -i^Le^'G(L^) 

we write the signal as it appears in the data as simply 

H{uj)^Gf{uj)hbar{i0) (19) 

with hiiar{^) given by the Fourier transform of Eq. |l^. 

At this point in the signal chain, in the interest of limiting the bandwidth of sample data required, the signal is 
mixed with a reference (whose frequency is chosen to be between the two detection mode frequencies) and low-pass 
filtered using a commercial lock-in amplifier |Q. The lock- in provides both in-phase {x) and quadrature {y) outputs, 
so both amplitude and phase information is available on a bandwidth of 125 Hz containing the resonant modes of the 
detector. In software, the in-phase and quadrature components were combined to form a complex data stream 

z{t) = x{t) + J y{t) 

which was then demodulted to a 1 Hz bandwidth around each of the resonant modes We describe the complete 
demodulation from 1 kHz to 1 Hz, including both the hardware and software lockins, as mixing the signal with a 
single reference exp [j{uj.t + 4>r)]j where is the reference frequency and (p^ is an unknown reference phase (this 
phase can in fact be measured, but for the data set in question was not). Returning to the time domain for clarity 
the mixed signal is 

H{t) = [Gf(t) * hbar{t)] exp []{uJrt + 0r)] (20) 

where indicates a convolution. Using the expression from Eq. ^ we get 

H(t) ^Gf{t)-k (i/ic (1 + cos^ i) [cos2V^o^^ (i) + sin2^/^o^x (*)] exp[-j(4>o - (j^r + {^o - ^r)t)] (21) 

I /ic (2cosi) [cos2V'oi^x W - sin 20^0 i^;(0] exp[-j($o - 0r + (wo - ^r)^)]) • 

Since LOr is on the order of 1 kHz, the effect of the low-pass filtering is to remove terms from H{t) which contain the 
sum frequencies {tOr + t^o)- Now returning to the frequency domain, the signal, after demodulation, can be written 

H{u) = i/ic(l + cos^i) exp[-j($o-0r)]G/(w)[cos2V'oi^+(w')+sin2-0oi^x(t^')] (22) 

- (2cosz) exp[-j($o - 0r)] G/(w) [cos 2-0o ^x (^^ ) - sin 2-00 F+(w )] 

where uj = ± {loq — i^r) is the (positive or negative) downconverted signal frequency. 
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VI. DETECTION CONSIDERATIONS 



Equation gives the form of the CW signal as it would appear in the Allegro data. We now ask: "Is this signal 
in the data?" Since the strength of the signal is small compared to the detector noise (otherwise we would see it on 
a spectrum analyzer!), some work needs to be done to answer this question. We will use the standard maximum 
likelihood method to guide the analysis, but in the end the question of whether a signal was detected or not will be 
answered experimentally. 

The likelihood function can be written 

with P{H\z) the probability that, given the observed data z{t), the signal H{t) is present, P{z\H) the probability 
that, given the signal H{t) is present, we observed the data z{t), and P{z\0) the probability that, given no signal 
present, we observe the data z{t). 
From H, if 

Z = Zl, Z2, ZN 

is the time-ordered sequence of detector output (in our case imaginary) and 

H = Hi,H2, ■■■■,Hn 
is the signal waveform at the same sample times, then 

1 1 

g,h — 

and 

1 / 1 ^"^ \ 

^(-'^)= [(4.)^^detiii^,.ii] -p -2 ^ ^;'^H^.-^.)(^-^'^r (25) 

where the Rgj^ is the inverse of the autocorrelation matrix. 



Substituting from above, the likelihood ratio is 

2 



A = exp(-- J2 Rghi-^gHr, - HgZ*h + HgH^)) (26) 



a, '1=0 

which, given the properties of the complex autocorrelation matrix can be written 

Af-l 1 , N-1 



A = exp hR R',rJ^9K " J E ^'ghH^K (27) 



2 

gji=0 J g,h=Q 



where 5R means "take the real part" . 

We define the following notation: Xk for the discrete Fourier transform of x(t), evaluated at uok = 2TTk/NAt, and 

Sn'u^ is the one-sided PSD. Using the relation pH] 



g,h=0 k=0 '^"'k 



E R~',x,xl = 2A/ E ^ (28) 



and taking the natural logarithm of the the likelihood function (to remove the exponential), we have 
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Substituting for H gives, 



lnA = 2A/5R(expH(<i>o + 0r)] 5] — ^ (30) 

X [^hc{l + cos^«) cos2V'oF^^ + (1 + cos^ i) sin2V^o^Xfc 

- ^ he (2 cos i) cos 2V'o ^'v + ^ (2 cos i) sin 2'(/'o ^"1 ] ) 
2 ''2 



fc=0 '-'"■fc 



X [^hlil + cos^ t)^ cos^2i;o\Gf,F+J' + ^hl{l+cos'i)' sin' 2^o\Gf,Fy^J 
+ \hl (2cosz)2 cos2 2MGf,FyA' + \hl (2cos»)2 sin^ 2MGf,F+,\'] 



where terms of the form 



fc=0 

sum to zero. 

To simphfy the notation of Eq. BOl we make the foUowing definitions: 



9+ 



^ zk G*f F* 

E /(I) (31) 



k=0 ^^k 
fc=0 '^"fc 



^^^^ ^n(^) 

fc=0 "fc 



-A/E 



fc=0 "fc 

ai = /ic(l + cos^ i) cos2?/'o 
^2 = hc{l + cos^ i) sin 2iI>q 
0-3, = ft.c(2 cosz) cos2?/'o 
04 = ^c(2 cos z) sin 2-0o • 

The four quantities gx j P+i Px } are famihar from signal processing theory. The first two, {qj^^q^}, are the 
outputs from applying independent optimal filters in the frequency domain for the two polarizations of the gravity 
wave. The second pair, {p+, px} are the signal to noise ratios (energy) for the two polarizations. These four terms are 
completely specified, up to the unknown signal frequency. The four a^'s are the "amplitudes" containing the desired 
astrophysical information. 

Substituting these into the likelihood function, 

In A = ai5R{exp[j($o - <l>r)] q+ } + a23? { exp[j($n - 0r)] Qx } 

- as^ij exp[j($o - (l)r)] qx} + a4^{j exp[j($o - (l^r)] q+ } 

- ^(a? P+ + (4px +alpy +alp+). 

We maximize the likelihood function over the 's which yields the following four expressions: 

ai = 2^^ cos($o - 0r + 4>+) (32) 
P+ 
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a2 =2l^cos($O"9!'r + 0x) 

03 = -2— sin($o - 0r + ) 

Px 

04 = 2— sin(i>o + (j)r- <^+) 

where the filtered ouputs have been expressed as 

9(+,x) = X) I exp(j X)) 

with (/)(+, X) = arg[g(+_x)]- Finally, 

^ = 4(%^ + = [(1 + cos2 z)2 + 4cos2 ^] = {hi + /i^ ) (33) 

which, to factors, is the energy of the gravitational wave. Note that this depends on the unknown signal frequency 
/o . As described below, we will assume values for the gravitational wave frequency, and use that to calculate Eq. |3|. 
The end result of the analysis is a "spectrum" of energy at a given signal frequency. We shall report this as a "strain 
amplitude" , given by the square root of Eq. 



hs{h) = ^h\{h) + hl{h). 



VII. EXPERIMENTAL CONSIDERATIONS 



Calculation of Eq. 33 involved a number of steps. First, the archived data was read off tape and narrowbanded to a 
1 Hz bandwidth around each of the resonant modes |9| where ALLEGRO is most sensitive. This was done to reduce 



the computational requirements. Next, the mode amplitudes were "cleaned" of large, transient event s (sect ion VII A) 
Finally, the cleaned data was discrete Fourier transformed and optimal filters were applied (section VII E|) . 



A. data selection 



Even with ALLEGRO'S high duty cycle, there were periods of missing or unusable data. Data losses came in 
basically three flavors. (1) Transient electronic effects which lasted on the order of a second. (2) Longer periods when 
the detector was undergoing some form of maintenance. (3) The clock losing phase lock to WWVB. The transient 
disturbances were the most frequent, occurring at a rate of one or two per day. They usually involved a sudden change 
in the flux threading the SQUID loop (hence the name "flux jumps"). The majority of the flux jumps occurred when 
the dc level of the SQUID reached a pre-determined maximum (5 volts). The electronics controlling the SQUID then 
reset the dc voltage to zero, causing a short and violent jump in the in-phase and quadrature channels of the data, as 
shown in Fig. |^. 

Frequently electronic interference reaching the SQUID caused flux jumps. In the past when data tapes were erased 
(the degaussing takes place in a separate room from the main experiment), the end of the degaussing cycle produced 
a noise spike which traveled through the wiring in the wall, through the computers, and from there to the SQUID. 
Once recognized, an isolating transformer was placed between the degausser and the wall socket, fixing the problem. 
Another common type of transient signal is a "spike" . These look similar to flux jumps in the data and there is some 
suspicion that they are in fact flux jumps, but essentially they are of unknown origin. All of these noisy periods were 
short enough so that the affected data could be removed and the resulting gap interpolated across. This was done 
using the MATLAB interpl routine. Interpolation was performed on the resampled data. It was usual to interpolate 
across 1-5 seconds of data, using the 10-20 seconds of data before and after the gap for the interpolation template. 

For longer sections of unusable data or for periods of missing data, the analysis was stopped and restarted after 
the disturbance. By "restarted" we mean that the accumulation of data was stopped, the accumulated data purged, 
and the analysis started up again on the data immediately following the disturbance. The most common cause of 
data loss was transferring liquid helium into the dewar which removed a couple of hours of data every week. Another 
cause of long stretches of unusable data were large excitations of the resonant modes due to earthquakes around the 
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globe. Earthquakes were identified by a unique signature in the low frequency housekeeping channel. It was usual for 
an earthquake to produce multiple large excitations over a few tens of minutes, often resulting in saturation of the 
A/D's. Computer down time, calibration of the detector and other maintenance all caused gaps in the data, although 
infrequently. 

The final type of data loss was associated with the WWVB clock. Frequently when the weather between Baton 
Rouge and NIST at Boulder, Colorado was bad, the clock we used to control the sampling of the data lost phase 
lock to the WWVB radio signal. When this happened, the clock's internal oscillator "freewheeled" with the result 
that the time between samples was no longer consistently 8 ms. Deviations in the sampling rate from 8 ms were 
called "timing jumps" . A jump in the time between when samples were taken produced a corresponding jump in 
the phase of modes and the calibrator signal as shown in Fig. ^ The most frequent jumps were on the order of 1-2 
ms, producing a phase jump in a sinusoidal signal at the mode frequencies of approximately 1/100 of a cycle. This 
was considered an acceptable level of uncertainty in the ability to track the phase of a potential gravity wave signal. 
These small jumps were noted but ignored. Larger jumps produced correspondingly larger jumps in phase and were 
considered unacceptable. When they occurred, the analysis was stopped at the timing jump and restarted again after 
the glitch. The frequency and size of the timing jumps were highly variable. During the winter of 1994 the smaller 
jumps occurred almost once per day while the larger jumps rarely happened. By the spring of that year, the trend 
was reversed and much data was lost due to the inconsistency of the clock. 



B. the filters 



The ability of our analysis filters to match the phase of the gravity wave determined the amount of continuous 
data which could be analyzed at one time. We refer to this as a "record" . Ideally, the length of a record would be 
set by the available computational facilities. This was not the case. Instead, the record length was limited by our 
electronics. The function generator which supplied the reference frequency to the lock-in detector had some drift 
associated with it. This drift limits the ability of the filter to match the phase of the gravity wave signal in the data. 
Careful measurements of the drift led to the conservative conclusion that roughly 28 hours (10^ sec) was the longest 
period of time for which we could expect the filter to remain in phase with the signal The less conservative 
conclusion was roughly three times longer, but we prefer to be cautious with the analysis. The reference signal to the 
lock-in has since been phase locked to GPS, greatly improving the phase stability. Given the above considerations, 
there were 34 records available from days 1-94 of 1994, for a total of 944.44 hours of data. 

Performing a discrete Fourier transform on 10^ sec worth of data results in a frequency resolution of the search 
of 10^^ Hz. The optimal filters also depend on the signal frequency. With no known source available, we assume a 
there is a potential signal every 10 /xHz in the ranges 896.30-897.30 Hz and 919.76-920.76 Hz. Each assumed signal 
frequency was used in turn to create the i<+ and components of the signal template. 

The other two components of the filters, the bar transfer function and the power spectral density, depend on the 
resonant frequency and damping time of the coupled bar-transducer system. As both quantities experienced slow 
drifts due to, for example, temperature changes in the dewar, it was necessary to measure them for each record (see 
Fig. H and Fig. ^ respectively). For each mode, a low variance power spectral density was formed from the cleaned 
data. A Lorenzian curve was then fit to each PSD using the MATLAB curvefit routine, which finds the best fit to a 
function in the least-squares sense. The Lorentzian was characterized by four parameters: the white-noise level (S'o), 
the peak height (Si), the frequency of the resonance (w±), and the width (or the damping time - t±). The resulting 
fitted parameters for each record were then stored on disk to be retrieved as necessary. The power spectral density 
component of the optimal filter was calculated for a record by 

5n(^,) =So+ ^-^^^[{u^l - u^ir + u?JtI\-^ . (34) 

Similarly, the bar transfer function was calculated from the functional 

G(a>fe) = fc± 5^ \ujl+2UJklT± - (35) 

where the parameters are those given above. The constant k± sets the calibration for each mode ||l7|| . 
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C. computational reduction 



Even though this search was directed towards only a small section of the sky, it was still cpu intensive on the 
available computing facilities. Applying 2 x 10^ different filters, where each new set of filter coefficients required a 
Fourier transform on a time sequence of 10^ elements, was prohibitively time consuming. 

The signal frequency enters into calculation of the filter coefficients in two places; through the carrier wave and in 
the Doppler shift. However, as the phase modulation is slowly varying with respect to the carrier wave, it was not 
necessary to re-calculate the phase correction due to the Doppler shift for each assumed signal frequency. Instead, the 
Doppler shift was calculated at specific choices of the signal frequency, which were then used to approximate signal 
templates over a range of frequencies. 

The maximum fractional frequency shift of a signal from the two sources over the course of a year is 

\Sf{fo)\ _ f 5 X 10-5 47 Tuc 

1 10^'* galactic center 

at the signal frequencies of interest. To approximate all phase shifts in a range of signal frequencies /q ± A/q with a 
calculation at a single signal frequency, then the error made must be less than the frequency resolution of the search. 
Defining the error as 

Sfifo ± A/o) - Sfifo) = 5 X lO-^A/o 

we have the maximum range of signal frequencies over which the approximation is valid is given by 

I < / 0.2 Hz, 47 Tuc 
I JO I _ <i Q 2^ Ylz, galactic center 

Using this information, we calculated independent F+, Fx components of the signal template at the frequencies 896.45 
Hz, 896.80 Hz, 897.15 Hz, 919.91 Hz, 920.26 Hz and 920.61 Hz for the 47 Tuc analysis. For the galactic center analysis, 
signal templates were calculated every 0.2 Hz in the range 896.40 Hz - 897.20 Hz for the minus mode, and 919.86 Hz 
- 920.66 Hz for the plus mode. 



VIII. RESULTS 



The real and imaginary components of the data stream, in the absence of a signal, are zero mean Gaussian 
distributed. Neither the Fourier transform or the filtering changes the underlying distribution and therefore the 
real and imaginary parts of the filtered outputs 9(+,x) are also individually zero mean Gaussian distributed. Given 
this, calculation of Eq. ^ at a particular choice of the signal frequency, taken from a particular data record, results 
in a sample drawn from a chi-squared distribution with four degrees of freedom |l8|. If the signal is absent, this is 
given by 

P{E) = exp {^E/2<j^) . (36) 

where E = hUfo). 

The total number of available data records provide an ensemble of experiments each making 2 x 10^ measures of E 
at the assumed values for /q. Since the signal is expected to be long-lived, much longer than the operational time of 
the detector, we improve our estimate of the energy by taking the ensemble average of the individual measures: 

1 ^ 

i=l 

where iV = 34 is the number of available data records. We calculate the strain amplitude /is(/o) by taking the square 
root of E. As each ensemble is generated by the sum of four terms (Eq. and there are 34 data records averaged 
together, from the central limit theorm we can accurately describe the distribution of /is(/o) as gaussian. 

The resulting "spectra" of hsifo) vs. fo are shown in Fig. ^ for 47 Tuc and Fig. ^ for the galactic center analysis. 
It is important to remember that the abscissa is not a Fourier frequency but rather the assumed frequency of the 
gravitational wave. That the plus mode shows a lower strain value than the minus mode reflects the fact that the 
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plus mode has a higher quahty factor and is therefore more sensitive. That the results of the galactic center analysis 
seem somewhat "shifted" with respect to the 47 Tuc results is due to the effects of the Doppler shifts. For a signal 
from 47 Tuc during the period over which data was analyzed, a signal at a particular frequency, for example 920 Hz, 
experienced periods of both red shift and blue shift. On average, therefore, the signal was experiencing detector noise 
at roughly 920 Hz. For the galactic center analysis, during this period an incoming signal was experiencing maximum 
blue-shift, and was never red-shifted. Therefore, a signal at 920 Hz would be experiencing detector noise at a slightly 
higher frequency, resulting in a graph that seems "a bit skewed to the right" with respect to the 47 Tuc analysis. 

The question still remains if there is a real signal at a particular frequency. It is unlikely that the detector is 
dominated by CW radiation. We expect the number of real signals which would be extracted by our analysis to be 
small with respect to the number of possible signals (i.e. the number of frequencies at which we assume there to be 
a signal). We can therefore answer the question experimentally. We normalize ^.^(/o) at each /p to its mean value. 
This new variable has a distribution independent of the detector noise at each signal frequency. 

The histograms of the resulting normalized spectra are shown in Fig. ^ for 47 Tuc and Fig. |lO| for the galactic center 
analysis. Given our reasonable assumption that the detector is not wave-dominated, the majority of the normalized 
measures form an experimental estimation of the parent distribution for the detector noise from which each individual 
measure is drawn. If the strain amplitude at a particular signal frequency lies outside this distribution, it is identified 
as a candidate for a real signal. 

IX. DISCUSSION 

Neither Fig. || or Fig. |l^ shows any evidence for "outliers" which would suggest the possibility of having observed 
gravitational radiation from pulsar spin-down. The existence of burst-like non-gaussian noise in resonant bar data 
is well known, and this noise must be dealt with in any search for burst-like signals of gravitational waves. The 
encouraging result of this analysis is that for CW signals there is not a serious problem with non-gaussian noise in the 
resonant bar data. Thus, at least at the level of the current work, all the noise sources are well understood and the 
overall sensitivity may be improved in the future simply by inceasing the length of the data record. We note that care 
was taken to use only data when the detector was operating well and procedures were implemeted to smooth out other 
non-stationary effects. Astrophysically, at this level of sensitivity, the actual spin-down rate of a pulsar due to the 
energy lost as gravitationl radiation would be large enough that the signal and its filter template would go out of phase 
by half a cycle after of order 5000 s, much less than the length of one record. This effect can be taken into account 
by including one or more spin-down parameters in the signal template, at the cost of significantly increasing the 
computational time of the search. Since there was no known source for ALLEGRO and somewhat limited computing 
resources available, we have essentially made a demonstration of the capabilities of resonant detectors for this type of 
search. 

However, it is important to note that even with these restrictions, the analysis has reached a level of sensitivity which 
is astrophysically interesting. If a source were to become known with the correct frequency, and, as our sensitivity 
was limited by hardware which has been substantially improved, it seems likely that such a signal could be detected 
and new astrophysics learned. 

X. CONCLUSION 

We have searched for sources of continuous gravitational radiation from the globular cluster 47 Tuc and the galactic 
center at signal frequencies near 1 kHz using data taken by the ALLEGRO gravitational radiation detector in the 
first three months of 1994. No candidate signals were found, and a constraint of 8 x 10"^"' was put on gravitational 
radiation emitted from pulsar spin-down at both locations. 
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FIG. 1. Frequency shift of a signal from 47 Tuc arriving at the ALLEGRO detector on January 1-2, 1994. 

FIG. 2. The reception patterns for the plus polarization (solid line) and the cross polarization (dot-dashed line) of a signal 
from 47 Tuc. 



FIG. 3. A SQUID reset in the in-phase and quadrature channels of the data stream. 

FIG. 4. A 1 ms jump in the sampling time as it appears in the continuous system test signal (CST-a sinusoidal driving force 
applied to the bar at 865 Hz) [Q. 

FIG. 5. The measured resonant frequencies for each mode and each data record. The freqeuncy is given as a shift relative 
to the measured frequency for the first data record. 



FIG. 6. The measured damping times for each mode and each record. 

FIG. 7. The strain amplitude of a gravitational wave from 47 Tuc at each assumed signal frequency near the plus and minus 
resonant modes. The lower values of the strains of the plus modes relative to the minus mode are due to the higher mechanical 
Q of the plus mode. 

FIG. 8. The strain amplitude of a gravitational wave from the galactic center at each assumed signal frequency near the plus 
and minus resonant modes. 

FIG. 9. Histograms of the normalized "spectrum" for signals from 47 Tuc near the minus and plus modes. The solid curve 
is a gaussian distribution with mean of one and variance derived from the normalized spectrum. 

FIG. 10. Histograms of the normalized "spectrum" for signals from the galactic center near the minus and plus modes. The 
solid curve is a gaussian distribution with mean of one and variance derived from the normalized spectrum. 
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TABLE I. Coordinates for the two source locations 



Source 
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